#include "linear_least_squares.h"

#include "math_util.h"
#include "matrix.h"

//
// Linear least squares method
//
double linear_least_squares(int degree, int num_elements, double* xi, double* yi, double* parameters)
{
	double** A = allocate_matrix(degree + 1, degree + 1);

	int max_power = degree * 2;
	for (int n = 0; n <= degree; n++)
	{
		for (int m = 0; m <= degree; m++)
			A[n][m] = sum_power(num_elements, xi, max_power - m);
		max_power--;
	}

	double** B = allocate_matrix(degree + 1, 1);
	for (int n = degree; n >= 0; n--)
		B[degree - n][0] = sum_mul_power(num_elements, xi, yi, n, 1);

	double** X = allocate_matrix(degree + 1, 1);
	solve_cramer(A, degree + 1, degree + 1, B, X);
	for (int i = 0; i < degree + 1; i++)
		parameters[i] = X[i][0];

	free_matrix(X, degree + 1, 1);
	free_matrix(B, degree + 1, 1);
	free_matrix(A, degree + 1, degree + 1);

	double error = 0.0;
	for (int i = 0; i < num_elements; i++)
	{
		double comp_y = evaluate_polynomial(xi[i], parameters, degree + 1);
		error += (yi[i] - comp_y) * (yi[i] - comp_y);
	}

	return error / double(num_elements);
}